# imports
import numpy as np
from aicspylibczi import CziFile
import matplotlib.pyplot as plt
import glob
import scipy.optimize
Background and Introduction¶
For our project, we are using image analysis techniques from BIOL300 to obtain data from images of FRAP (fluorescence recovery after photobleaching) techniques over time. We can use the data along with techniques from BIOL301 to model the diffusion behavior and recovery curve of a specific FRAP dataset that we will choose.
For this project, we chose to use FRAP data from Dr. Nikki Farnsworth's paper Modulation of Gap Junction Coupling Within the Islet of Langerhans During the Development of Type 1 Diabetes. This paper examines the role of gap junctions, specifically Connexin 36 (Cx36), in Type 1 Diabetes (T1D). Therefore, the primary goal was to determine if changes in islet Cx36 gap junctions coupling and $Ca^{2+}$ signaling dynamics occur early in the progression of T1D, prior to changes in glucose tolerance. FRAP was used in this research paper so that one would be able to measure the rate at which fluorescently labeled molecules bound to Cx36 diffused between cells. This allowed for the comparison of FRAP recovery rates in islets from healthy versus diabetic mice. The FRAP results help visualize how altering gap junction coupling can affect insulin secretion and $Ca^{2+}$ signaling.
The intellectual merit of this project involves applying models of box-jumping plots and differential equations to the experimental technique of FRAP to graph and analyze FRAP-related data. FRAP can be a highly useful tool for examining movement in cells, including the diffusion of molecules, lateral and transverse movement of phospholipids in the cell membrane, and other processes. By modeling this graphically and mathematically, we can gain a better understanding of the data and the biological process represented by it.
The broader impact of this project is that FRAP analysis has the potential to make significant contributions to various fields including biological research, drug discovery, and environmental science. This project has the ability to provide insights into protein-protein interactions, cellular signaling, and enzyme activity. In addition, it can help increase our understanding of processes such as diffusion, transportation, and other cellular behaviors. Next, FRAP analysis can help to identify proteins involved in diseases by studying their mobility in cells. We can also use this information to see the impact of certain drugs on protein-protein interactions and cellular processes. Finally, FRAP analysis can be used to increase understanding on the impact of certain pollutants on cellular processes and organism health.
Defining Functions¶
extract_czi()¶
The purpose of the extract_czi() function is to process CZI image files and extract a sequence of image frames from them, each image frame representing a specific time point during the FRAP experiment. Then, these image frames will be organized into a list of image data by time.
def extract_czi(path):
'''Given a path to a czi file, reads data from the file into a list of images by time.'''
# read the czi file
czi = CziFile(path)
dims = czi.get_dims_shape()[0] # dictionary of dimensions (X, Y, and T are the ones we care about)
# get time dimension
time = dims['T']
# get the array of image data at each time point
images = []
for t in range(time[1]):
im = czi.read_image(T=t)
images.append(im[0])
# remove extra dimensions of each image (squeeze)
for i in range(len(images)):
images[i] = np.squeeze(images[i])
return images
# see Image Analysis subsection for testing
graph_subplots()¶
The purpose of the graph_subplots() is to visualize a list of images as a grid of subplots. This function will enable us to see the time series data.
def graph_subplots(data, num_rows, num_columns, title=True):
'''Takes image data and graphs it as subplots.'''
# set axes
fig, axs = plt.subplots(num_rows, num_columns, figsize=(15,15))
plt.tight_layout()
# graph each image
for i, ax in enumerate(axs.flat):
if i < len(data):
ax.imshow(data[i])
# label with time step if desired
if title:
ax.set_title(f"Time Step t={i+1}")
# see simulate diffusion test case for testing
slice_image()¶
The purpose of the slice_image() function is to help extract a rectangular portion from a given image.
def slice_image(image, x, y, width, height):
'''Slices and returns a rectangular portion of the given image.'''
# initialize list to store sliced image data
sliced = []
# add each pixel in the region to slice
for i in range(y, y+height):
new_row = []
for j in range(x, x+width):
new_row.append(image[i][j])
sliced.append(new_row)
return np.array(sliced)
# see Image Analysis subsection for testing
average_fluorescence()¶
The purpose of this function is to help quantify the overall fluorescence intensity in a given region of interest over time. The function does this by iterating through each image in the input list and then calculates the average pixel intensity for each row. The function then returns the list of average fluorescence values.
def average_fluorescence(images):
'''Given a list of images, returns a list of the average fluorescence for each.'''
# initialize list to store average fluorescences
avg_fluorescences = []
# iterate through each image
for im in images:
im_fluorescence = []
# store the average fluorescence of each row
for row in im:
im_fluorescence.append(np.mean(row))
# find the total average fluorescence of the image and add it to the list of average fluorescences
avg_fluorescences.append(np.mean(im_fluorescence))
return avg_fluorescences
# see Fluorescence of Unbleached vs Bleached Regions for testing
find_nearest_pixel()¶
The purpose of the find_nearest_pixel() is to help find the nearest fluorescence pixel in a specific direction from the current pixel.
def find_nearest_pixel(image, location, direction, fluorescable_map):
'''Given an image and the location of a pixel in that image, returns tuple (x,y) of the
location of the nearest fluorescable pixel within an 11x11 block in a given direction,
or -1 if a pixel cannot be found.'''
# extract x and y coordinates from location
x_loc = location[0]
y_loc = location[1]
# build the 11x11 block to check and a corresponding fluorescence map
# initialize lists to store the block and fluorescence map
block = []
block_map = []
# find indices of the edges of the block
left_edge = x_loc - 5
if left_edge < 0: # in case block hits the bounds of the image
left_edge = 0
right_edge = x_loc + 6
if right_edge > len(image[0]):
right_edge = len(image[0])
up_edge = y_loc - 5
if up_edge < 0:
up_edge = 0
down_edge = y_loc + 6
if down_edge > len(image):
down_edge = len(image)
# create the block by adding the pixels from the original image
y_count = 0
for y in range(up_edge, down_edge):
block_row = []
map_row = []
x_count = 0
for x in range(left_edge, right_edge):
# store the new location of the pixel we're searching from
if x == x_loc and y == y_loc:
block_x = x_count
block_y = y_count
# add the corresponding pixel from the image to the block and the map
block_row.append(image[y][x])
map_row.append(fluorescable_map[y][x])
x_count += 1
block.append(block_row)
block_map.append(map_row)
y_count += 1
# initialize variables to store the minimum directional distance, minimum total distance, and nearest pixel location
min_dist = np.inf
min_total = np.inf
nearest_p = -1
# find nearest pixel within block
for y, row in enumerate(block):
for x, p in enumerate(row):
# skip any pixels that aren't fluorescable based on the fluorescence map
if not block_map[y][x]:
continue
# calculate directional and total distances
left_dist = block_x - x
right_dist = -left_dist
up_dist = block_y - y
down_dist = -up_dist
total_dist = (left_dist**2 + up_dist**2)**0.5
# depending on selected direction, stores the location of the pixel if it's closer than the existing closest
if direction == "left":
if left_dist <= min_dist and left_dist > 0 and total_dist <= min_total:
nearest_p = x, y
min_dist = left_dist
min_total = total_dist
# catch left border
if location[0] == 0:
return -1
elif direction == "right":
if right_dist <= min_dist and right_dist > 0 and total_dist <= min_total:
nearest_p = x, y
min_dist = right_dist
min_total = total_dist
# catch right border
if location[0] == len(image[0]):
return -1
elif direction == "up":
if up_dist <= min_dist and up_dist > 0 and total_dist <= min_total:
nearest_p = x, y
min_dist = up_dist
min_total = total_dist
# catch top border
if location[1] == 0:
return -1
elif direction == "down":
if down_dist <= min_dist and down_dist > 0 and total_dist <= min_total:
nearest_p = x, y
min_dist = down_dist
min_total = total_dist
# catch bottom border
if location[1] == len(image):
return -1
# convert coordinates of nearest pixel to be relative to the image rather than the block
if nearest_p != -1:
x_dist = nearest_p[0] - block_x
y_dist = nearest_p[1] - block_y
nearest_x = x_loc + x_dist
nearest_y = y_loc + y_dist
nearest_p = nearest_x, nearest_y
return nearest_p
# test the block construction of find nearest pixel function
def test_check_block(x_loc, y_loc, image):
# build the 11x11 block to check and a corresponding fluorescence map
# initialize lists to store the block and fluorescence map
block = []
block_map = []
# find indices of the edges of the block
left_edge = x_loc - 5
if left_edge < 0: # in case block hits the bounds of the image
left_edge = 0
right_edge = x_loc + 6
if right_edge > len(image[0]):
right_edge = len(image[0])
up_edge = y_loc - 5
if up_edge < 0:
up_edge = 0
down_edge = y_loc + 6
if down_edge > len(image):
down_edge = len(image)
# create the block by adding the pixels from the original image
y_count = 0
for y in range(up_edge, down_edge):
block_row = []
# map_row = [] ignoring map for test case
x_count = 0
for x in range(left_edge, right_edge):
# store the new location of the pixel we're searching from
if x == x_loc and y == y_loc:
block_x = x_count
block_y = y_count
# add the corresponding pixel from the image to the block and the map
block_row.append(image[y][x])
# map_row.append(fluorescable_map[y][x])
x_count += 1
block.append(block_row)
# block_map.append(map_row)
y_count += 1
return block
image = [
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
]
x_loc = 6
y_loc = 7
test_check_block(x_loc, y_loc, image)
[[1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]]
# test find nearest pixel function
test_image = [
[0, 0, 0, 1, 1],
[0, 0, 0, 0, 0],
[0, 0, 1, 0, 0],
[0, 0, 0, 0, 0],
[1, 0, 0, 1, 0]
]
test_map = [
[0, 0, 0, 1, 1],
[0, 0, 0, 0, 0],
[0, 1, 1, 0, 0],
[0, 0, 0, 0, 0],
[1, 0, 0, 1, 0]
]
test_loc = (2, 0) # 0 in middle of first row
print(f"Left: {find_nearest_pixel(test_image, test_loc, 'left', test_map)}") # expected (1, 2)
print(f"Right: {find_nearest_pixel(test_image, test_loc, 'right', test_map)}") # expected (3, 0)
print(f"Up: {find_nearest_pixel(test_image, test_loc, 'up', test_map)}") # expected -1
print(f"Down: {find_nearest_pixel(test_image, test_loc, 'down', test_map)}") # expected (2, 2)
Left: (1, 2) Right: (3, 0) Up: -1 Down: (2, 2)
simulate_diffusion()¶
The purpose of simulate_diffusion() is to help simulate the process over time given an initial image and fluorescence map. This will help us utilize a probabilistic approach to model the movement of the fluorescent molecules between pixels.
def simulate_diffusion(image, fluorescable_map, delta_t, num_steps, k):
'''Given a starting image and map of fluorescable pixels, simulates diffusion using probability and returns
array of image data to display the results of the simulation.'''
# get height and width of image
height = len(image)
width = len(image[0])
# build and initialize probability array
p = np.zeros((num_steps, height, width))
# initial condition is the fluorescence of our starting image--more fluorescent means it's more likely to diffuse
for y in range(height):
for x in range(width):
p[0][y][x] = image[y][x]
# model diffusion over time
for t in range(num_steps-1):
# loop over each pixel
for y in range(height):
for x in range(width):
# don't diffuse if the pixel can't fluoresce
if fluorescable_map[y][x] == 0:
continue
pix_loc = (x, y)
# find the nearest pixel in each direction--these are the "boxes" fluorescence can diffuse to
left = find_nearest_pixel(image, pix_loc, "left", fluorescable_map)
right = find_nearest_pixel(image, pix_loc, "right", fluorescable_map)
up = find_nearest_pixel(image, pix_loc, "up", fluorescable_map)
down = find_nearest_pixel(image, pix_loc, "down", fluorescable_map)
# store each of the valid pixels
boxes = [left, right, up, down]
diffusable_boxes = [box for box in boxes if box != -1]
# initialize new probability
new_prob = 0
# add to the new probability for fluorescence diffusing into our current pixel
for box in diffusable_boxes:
new_prob += p[t, box[1], box[0]]
# remove from the new probability for fluorescence diffusing out of our current pixel
new_prob -= len(diffusable_boxes) * p[t, y, x]
# update probability matrix
p[t+1, y, x] = p[t, y, x] + k*delta_t*new_prob
return p
The following is a test case of our simulate_diffusion() function to show how an initial image will demonstrate modeled FRAP behavior over time given a fluorescence map. As seen below, the pixels with high intensity at first grow dimmer over time and the initially dark pixels grow brighter over time as the "fluorescent molecules" diffuse into them. This is the behavior we would expect to see for this model.
# test diffusion simulation
t_steps = 20
dt = 0.1
k_guess = 1
test_image = [
[0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0],
[1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0],
[1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1],
[0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1],
[1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0]
]
test_map = [
[0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0],
[1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0],
[1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1],
[0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1],
[1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1],
[0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0],
[1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1],
[1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1],
[0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0]
]
p = simulate_diffusion(test_image, test_map, dt, t_steps, k_guess)
graph_subplots(p, 5, 4)
CZI class¶
The purpose of the CZI class is to help with future applications of this project by allowing us to organize different attributes of a czi file into one object, making it easy to repeat steps taken in this project with other FRAP data. By using the set_bleached_dims() and set_unbleached_dims() functions allows one to define separate bounding boxes for bleached and unbleached regions within the chosen (FRAP czi) image. Then, the slice_image() function iterates through each image and slices it based on the defined region. This stores the unbleached, bleached, and combined data in separate lists.
class CZI:
def __init__(self, path):
'''Constructs CZI object using file path'''
# store critical information
self.path = path
self.images = extract_czi(self.path)
self.time_steps = len(self.images)
# initialize dimensions for bleached and unbleached slices
self.x_b = None
self.y_b = None
self.width_b = None
self.height_b = None
self.x_u = None
self.y_u = None
self.width_u = None
self.height_u = None
self.bleach_ind = None
# initialize lists for sliced regions
self.sliced_data = []
self.unbleached_data = []
self.bleached_data = []
def set_bleach_ind(self, bleach_ind):
'''Sets and stores the time index at which photobleaching occurred.'''
self.bleach_ind = bleach_ind
def set_bleached_dims(self, x, y, width, height):
'''Sets and stores the dimensions and location of the bleached region of the images.'''
self.x_b = x
self.y_b = y
self.width_b = width
self.height_b = height
def set_unbleached_dims(self, x, y, width, height):
'''Sets and stores the dimensions and location of the unbleached region of the images.'''
self.x_u = x
self.y_u = y
self.width_u = width
self.height_u = height
def slice_images(self):
'''Uses bleached and unbleached dimensions to slice the image data and store it in the object.'''
# loop over each image and slice it
for im in self.images:
unbleached_slice = slice_image(im, x_u, y_u, width_u, height_u)
bleached_slice = slice_image(im, x_b, y_b, width_b, height_b)
# store sliced data
self.unbleached_data.append(unbleached_slice)
self.bleached_data.append(bleached_slice)
self.sliced_data.append(np.concatenate((unbleached_slice, bleached_slice)))
def show_images(self):
'''Displays every image in the czi file.'''
to_graph = []
for image in self.images:
to_graph.append(image.data)
graph_subplots(to_graph, 5, 5)
def __str__(self):
'''Establishes string to display when object is printed.'''
return self.path
Image Analysis¶
To conduct image analysis, we began by inputting Dr. Farnsworth's FRAP czi data from her study Modulation of Gap Junction Coupling Within the Islet of Langerhans During the Development of Type 1 Diabetes. Once the data was input, we then extracted the image data from every czi in the folder. To make sure our code was working properly, we tested the first czi file in the folder, RIP_0.01_FRAP1.czi, where RIP stands for Rat Insulin Promoter. It is important to note that the following code can easily be applied to any of the images in the folder, it would just require adjusting the slices and indices.
# get files
path_list = glob.glob("FRAP_data/02102017/*.czi")
path_list
['FRAP_data/02102017\\RIP_0.01_FRAP1.czi', 'FRAP_data/02102017\\RIP_0.01_FRAP2.czi', 'FRAP_data/02102017\\RIP_0.01_FRAP3.czi', 'FRAP_data/02102017\\RIP_0.1_FRAP1.czi', 'FRAP_data/02102017\\RIP_0.1_FRAP2.czi', 'FRAP_data/02102017\\RIP_0.1_FRAP3.czi', 'FRAP_data/02102017\\RIP_0_FRAP1.czi', 'FRAP_data/02102017\\RIP_0_FRAP2.czi', 'FRAP_data/02102017\\RIP_0_FRAP3.czi', 'FRAP_data/02102017\\WT_0.01_FRAP1.czi', 'FRAP_data/02102017\\WT_0.1_FRAP1.czi', 'FRAP_data/02102017\\WT_0.1_FRAP2.czi', 'FRAP_data/02102017\\WT_0_FRAP1.czi', 'FRAP_data/02102017\\WT_0_FRAP2.czi', 'FRAP_data/02102017\\WT_0_FRAP3.czi']
# store each file as a czi object
czi_files = []
for path in path_list:
czi_files.append(CZI(path))
# select czi file to work with and display
czi = czi_files[0]
czi.show_images()
The above images are our dataset which captures the fluorescence intensity over time. Each time step in this czi file is 15 seconds apart.
The images will then be sliced into two regions: bleached (which will initially be dark after photobleaching) and unbleached (which will retain the fluorescence) so we can examine the behavior of the separate regions.
# define coordinates for unbleached region to slice
x_u = 98
y_u = 82
width_u = 307
height_u = 180
czi.set_unbleached_dims(x_u, y_u, width_u, height_u)
# define coordinates for bleached region to slice
x_b = 98
y_b = 262
width_b = 307
height_b = 215
czi.set_bleached_dims(x_b, y_b, width_b, height_b)
# slice the images and store
czi.slice_images()
unbleached_data = czi.unbleached_data
bleached_data = czi.bleached_data
sliced_data = czi.sliced_data
# display sliced images
graph_subplots(sliced_data, 5, 5)
The above sliced images are a more helpful way to show fluorescence recovery of Cx36, as they are zoomed in on the cell. After photobleaching, which occurs at time step t=4, the increase in fluorescence intensity in the bleached region as time passes is consistent with the expected recovery due to molecular diffusion. As well as this, the decrease in fluorescence intensity in the unbleached region corresponds to molecules diffusing out to restore the equilibrium. These images show that molecular diffusion is occurring as expected.
We can then zoom in on the image at t=4 to see the sliced bleached and unbleached regions separately, as well as the full sliced image, below:
# show unbleached, bleached, and total sliced regions immediately after photobleaching has occurred
czi.set_bleach_ind(3)
bleach_ind = czi.bleach_ind
slice_zooms = [unbleached_data[bleach_ind], bleached_data[bleach_ind], sliced_data[bleach_ind]]
graph_subplots(slice_zooms, 1, 3, title=False)
Fluorescence of Unbleached vs Bleached Regions¶
In the below graph, the average fluorescence is calculated for both the unbleached and bleached regions for each time step. From the comparison of the average fluorescence, we can see that the fluorescence of the unbleached region decreases slightly over time as molecules diffuse out. In the bleached region, we can see that the fluorescence drops suddenly at the time step of the bleaching, then recovers as time goes on. This graph aligns with the expected behavior of FRAP experiments, which helps to confirm the diffusion of fluorescent molecules from the unbleached region to the bleached region over time.
# calculate average fluorescences of unbleached and bleached regions
unbleached_fluorescences = average_fluorescence(unbleached_data)
bleached_fluorescences = average_fluorescence(bleached_data)
# get time values
t = np.arange(czi.time_steps)
# graph average fluorescence over time for each region
plt.plot(t, bleached_fluorescences, '.', label="Bleached Region")
plt.plot(t, unbleached_fluorescences, '.', label="Unbleached Region")
plt.xlabel("Time")
plt.ylabel("Average Fluorescence")
plt.legend(loc="right")
plt.show()
Simulated Diffusion¶
One of the ways we modeled FRAP was with a probability-based simulation. The concept is derived from the box-jumping master equation model, where we create a matrix of the probability a particle is in a given “box” and update that matrix over a set number of time steps. We expanded this model from a linear one-dimensional model to a two-dimensional model so we could use it to model the diffusion in this two-dimensional cell.
We took the first image after photobleaching occurred and thresholded it to have binary fluorescence values (0 or 1), then used that as the start image for our simulation. The simulation then, for each pixel in the image, finds the nearest “fluorescable” pixel above, below, to the left, and to the right of that pixel. These pixels are determined to be fluorescable based on a fluorescence map, which is another thresholded image from the start of the experiment. The fluorescence of each of these pixels is directly related to how likely it is to diffuse, so the simulation uses probability to adjust the fluorescence of each. It does this for every time step, then we can examine the grid of images it generates to see the simulation of diffusion over time.
# create a fluorescability map by thresholding the first image in the file
map_base = sliced_data[0]
thresh = 4000
fluorescable_map = map_base > thresh
plt.imshow(fluorescable_map);
# create an image to start the diffusion simulation by thresholding the image immediately after photobleaching has occurred
start = sliced_data[bleach_ind]
thresh = 4000 # same threshold as map for consistency
sim_start = start > thresh
plt.imshow(sim_start);
# set parameters
t_steps = czi.time_steps - bleach_ind # run simulation for as many time steps as the czi has after photobleaching
dt = 0.1
k_guess = 1
# run diffusion simulation and display results
p = simulate_diffusion(sim_start, fluorescable_map, dt, t_steps, k_guess)
graph_subplots(p, 5, 4)
# compare beginning and end of diffusion simulation side by side
to_compare = (p[0], p[-1])
graph_subplots(to_compare, 1, 2, title=False)
In the images above, which are the beginning and end of the simulation respectively, we can see that fluorescence has spread out through the simulated cell to start filling gaps, as would be expected in the diffusion process. These can then be compared to the beginning and end of diffusion in the original image data, which is shown below:
# show beginning and end of sliced data to compare
to_compare = (sliced_data[0], sliced_data[-1])
graph_subplots(to_compare, 1, 2, title=False)
Differential Equations¶
The second way we modeled FRAP was with a plot of the integrated differential equation $\frac{dC}{dt} = k_{off}(C_{tot} - C) - k_{on}C$, which describes FRAP behavior. The goal was to use this as our predicted curve and compare it to our actual curve. $C$ will be our fluorescence values over time, $C_{tot}$ will be the total concentration of fluorescent molecules, $k_{off}$ will be the dissociation rate constant, and $k_{on}$ will be the association rate constant.
First, we have created the FRAP_rhs function to define the right-hand side of the differential equation that we will be integrating to find our predicted FRAP curve.
def FRAP_rhs(C, t, C_tot, k_off, k_on):
'''Describes right-hand side of FRAP differential equation.'''
# compute dC/dt
dC_dt = k_off*(C_tot - C) - k_on*C
# Return the result as a NumPy array
return np.array(dC_dt)
After defining our function, we can plot our predicted and actual data over time on the same graph. We will the average_fluorescence function to get out actual data. This will help us visualize the behavior of our actual FRAP data versus what we would expect it to do. We first defined the initial values and parameters of the differential equation we are going to integrate. For now, some of these will be guesses that we can later find optimized values for. Then, using the scipy.integrate.odeint function, we integrated our DE to create our predicted curve (C_guess). The predicted curve and the actual data are plotted against time on the same plot below. We can see visually that the predicted curve fits our data fairly well.
# integrate diff eq and graph with guesses
# get data (average fluorescences after bleaching) and time points
data = np.array(average_fluorescence(czi.images[bleach_ind:]))
t = np.arange(len(data))
# initial condition guess
C_0_guess = data[0]
# parameter guesses
C_tot_guess = data[0]
k_off_guess = 0.15
k_on_guess = -0.008
# package parameters into a tuple
guess_args = (C_tot_guess, k_off_guess, k_on_guess)
# integrate diff eq
guess_fit = scipy.integrate.odeint(FRAP_rhs, C_0_guess, t, args=guess_args)
# compare actual data to predicted curve found by numerical integration
plt.plot(t, data, '.', label="Data")
plt.plot(t, guess_fit, label="Guess Fit")
plt.xlabel("Time")
plt.ylabel("Fluorescence")
plt.legend()
plt.show()
After doing this with "best guess" values for some of our parameters, we then tried to find optimized values. Below we defined a function to return the chi-squared value for our provided parameters and data. We then used this function to perform scipy.optimize.minimize in order to find the optimal values we should use based on our guesses.
def chi_squared_FRAP(params, data, t):
'''Returns the chi squared value for the parameters and data provided.'''
# unpack the parameters (guesses)
C_tot = params[0]
k_off = params[1]
k_on = params[2]
C_0 = params[3]
arguments = (C_tot, k_off, k_on)
# compute the fit by integrating diff eq
fit = scipy.integrate.odeint(FRAP_rhs, C_0, t, args=arguments)
# compute the chi-squared
vals = (data - fit)**2 / fit
return np.sum(vals)
# optimize values by minimizing chi squared
optimal_vals = scipy.optimize.minimize(chi_squared_FRAP, [C_tot_guess, k_off_guess, k_on_guess, C_0_guess], args=(data, t),
bounds=[ (2000, 2500), (-1, 1), (-1, 1), (2000, 2500) ])
# extract optimal values
best_C_tot = optimal_vals.x[0]
best_k_off = optimal_vals.x[1]
best_k_on = optimal_vals.x[2]
best_C_0 = optimal_vals.x[3]
print(f'Best C tot: {best_C_tot}, Best k off: {best_k_off}, Best k on: {best_k_on}, Best C_0: {best_C_0}')
Best C tot: 2006.3229817087017, Best k off: 0.9999902430202654, Best k on: -0.0712297307784808, Best C_0: 2160.1942355373376
After finding the optimal values, we plugged them back into our DE and once again integrated to find the curve of scipy fit. We plotted this curve over time alongside our actual data and our initial guess fit curve. As seen below, this plot shows how the scipy fit compares to the guess fit. The results were unexpected, as something went wrong with the the scipy.optimize.minimize function and resulted in a flat line. Ideally, the scipy fit would have resulted in a curve that more closely matches the actual data than the guess fit, but that did not happen in this case. So instead, we used np.polyfit to calculate a polynomial function that should fit the actual data. This curve was also graphed on the plot below and does in fact match the actual data more closely. This method demonstrates how a FRAP curve can be predicted using an optimization integration of this particular FRAP differential equation model.
# calculate best fit by integrating with optimized values
best_args = (best_C_tot, best_k_off, best_k_on)
best_fit = scipy.integrate.odeint(FRAP_rhs, best_C_0, t, args=best_args)
# calculate polynomial fit
coeffs = np.polyfit(t, data, 2)
poly_func = np.poly1d(coeffs)
poly_fit = poly_func(t)
# plot data and fit curves
plt.plot(t, data, '.', label="Data")
plt.plot(t, guess_fit, label="Guess Fit")
plt.plot(t, best_fit, label="Scipy Fit")
plt.plot(t, poly_fit, label="Poly Fit")
# label and show graph
plt.xlabel('Time')
plt.ylabel('Fluorescence')
plt.legend()
plt.show()
Results and Discussion¶
We were able to successfully utilize functions in order to extract a series of images over time from a czi file as well as a list of average fluorescences. These functions can be used for any czi file containing FRAP experimental data and they worked the way we expected them to. The images and plotted data we got from the czi file of FRAP data that we used served as the basis for our models.
Next, we defined functions to create our first model: a box-hopping plot that could emulate the fluorescence recovery behavior of a cell after being photobleached. This model was able to successfully show how initially intense pixels from the unbleached region would diffuse into the bleached region over time, causing those pixels to become brighter. This model requires an image of an actual cell that has just undergone photobleaching and a given fluorescence map to work, but the subsequent modeled behavior over time emulates the actual FRAP data well. A model like this could provide useful information about how long a specific cell is predicted to take before fluorescence is recovered. This could be determined depending on how many time steps it takes for the image to display an equal amount of fluorescent pixels in all regions.
The second model we defined relied on functions for a differential equation of FRAP behavior, an integration of this DE, and an optimized polynomial fit. We first demonstrated that an integration of the differential equation would produce a curve that actually fit our data. We used guess values at first for the coefficients and this resulted in a fairly good fit of the actual data. We then attempted to use scipy.optimize to make a better fit by finding optimal values for the coefficients; however, this function did not produce the results we expected and did not work correctly. Since this method did not work to create a better fit, we instead used a function to create an optimized polynomial curve. This fit our actual data more closely than the guess fit we initially created, making it a good method to use for plotting FRAP behavior over time. This model could be useful for predicting what the recovery curve of a FRAP experiment under certain conditions would be. These conditions could be used to define the coefficients of the DE, and based on these values, this model could predict the optimal curve of fluorescent molecules recovering throughout a cell over time.
Overall, the models we created behaved mostly as expected and were ultimately able to match the actual data we were using fairly well. While there were a few problems we encountered, we were able to work around them to successfully create methods of analyzing and modeling FRAP data from a file.
Next Steps¶
If we were to explore further methods of analyzing FRAP data, we could start by comparing the FRAP curves of cells under different conditions. For instance, the dataset we used to for this project contained multiple files of FRAP data that were measured under different conditions. The purpose was the track the diffusion of Cx36 between the cells of mice that were genetically edited to display different diabetic conditions. We could use the methods of modeling that we developed to predict the behavior of Cx36 diffusion under each set of conditions. When comparing this to the actual data, the results should match closely to support whatever effect is observed.
We could also expand our analysis methods to use them on data from a single cell in which multiple FRAP experiments were performed. By labeling different types of molecules in areas of interest throughout the cell with different fluorescent colors, multiple FRAP experiments could be run at the same time. We would need to observe the diffusion behavior of multiple different molecules over time after photobleaching, likely in different areas of the same cell. We could isolate each molecule of interest by color and then use the same image slicing methods from this analysis to isolate the area of interest in the cell where the molecule is located. Developing ways to analyze complex data like this could allow us to model the diffusion behavior of multiple molecules at once. This is promising for research into how the behavior of one molecule would affect the diffusion of the other.